# -*- coding: utf-8 -*-
"""
    功能：基于遗传算法确定最优违约触发点
    版本：v1.0
"""

from scipy import stats
from scipy.optimize import fsolve
import numpy as np
import pandas as pd
import geatpy as ga

# 公司总样本数 / 2
N = 64
# ST 和非 ST 公司数均为 N / 2
num_com = int((N / 2))
# ST 和 ST*公司数均为
count = int((num_com / 2))

global counter
counter = 0

# 读取ST上市公司数据
ST_data_xls_1 = pd.read_excel("ST公司数据.xls")
ST_data_xls_2 = pd.read_excel("ST公司数据负债.xls")
# 读取ST*上市公司数据
ST_Plus_data_xls_1 = pd.read_excel("ST+公司数据.xls")
ST_Plus_data_xls_2 = pd.read_excel("ST+公司数据负债.xls")
# 读取非ST上市公司数据
NST_data_xls_1 = pd.read_excel("非ST公司数据.xls")
NST_data_xls_2 = pd.read_excel("非ST公司数据负债.xls")

# 读取ST上市公司数据
ST_data_1 = ST_data_xls_1.ix[: count,
            ['收盘价_Clpr', '流通股_Trdshr', '已上市流通股_Lsttrdshr', '年收益率_Yrret', '年无风险收益率_Yrrfret', '每股净资产(元/股)_NAPS']].values
ST_data_2 = ST_data_xls_2.ix[: count,
            ['流动负债合计(元)_TotCurLia', '非流动负债合计(元)_TotNCurLia', '负债合计(元)_TotLiab']].values

# 读取ST*上市公司数据
ST_Plus_data_1 = ST_Plus_data_xls_1.ix[: count,
                 ['收盘价_Clpr', '流通股_Trdshr', '已上市流通股_Lsttrdshr', '年收益率_Yrret', '年无风险收益率_Yrrfret',
                  '每股净资产(元/股)_NAPS']].values
ST_Plus_data_2 = ST_Plus_data_xls_2.ix[: count,
                 ['流动负债合计(元)_TotCurLia', '非流动负债合计(元)_TotNCurLia', '负债合计(元)_TotLiab']].values

# 读取非ST上市公司数据
NST_data_1 = NST_data_xls_1.ix[: num_com,
             ['收盘价_Clpr', '流通股_Trdshr', '已上市流通股_Lsttrdshr', '年收益率_Yrret', '年无风险收益率_Yrrfret', '每股净资产(元/股)_NAPS']].values
NST_data_2 = NST_data_xls_2.ix[: num_com,
             ['流动负债合计(元)_TotCurLia', '非流动负债合计(元)_TotNCurLia', '负债合计(元)_TotLiab']].values

# ST公司数据
# 收盘价_Clpr
ST_Clpr_list = ST_data_1[:, 0]
# 流通股_Trdshr
ST_Trdshr_list = ST_data_1[:, 1]
# 已上市流通股_Lsttrdshr
ST_Lsttrdshr_list = ST_data_1[:, 2]
# 每股净资产
ST_NAPS_list = ST_data_1[:, 5]

# 无风险列表
ST_r_list = ST_data_1[:, 4]
# 股权收益率波动率列表
ST_PriceTheta_list = ST_data_1[:, 3]
# 股权市场价值列表
ST_E_list = ST_Clpr_list * ST_Lsttrdshr_list + (ST_Trdshr_list - ST_Lsttrdshr_list) * ST_NAPS_list
# Short-term Debt (in 10,000)
ST_SD_list = ST_data_2[:, 0]
# print(ST_SD_list)
# Long-term Debt (in 10,000)
ST_LD_list = ST_data_2[:, 1]
# print(ST_LD_list)
# 负债合计
ST_D_list = ST_data_2[:, 2]

# ST*公司数据
# 收盘价_Clpr
ST_Plus_Clpr_list = ST_Plus_data_1[:, 0]
# 流通股_Trdshr
ST_Plus_Trdshr_list = ST_Plus_data_1[:, 1]
# 已上市流通股_Lsttrdshr
ST_Plus_Lsttrdshr_list = ST_Plus_data_1[:, 2]
# 每股净资产
ST_Plus_NAPS_list = ST_Plus_data_1[:, 5]

# 无风险列表
ST_Plus_r_list = ST_Plus_data_1[:, 4]
# 股权收益率波动率列表
ST_Plus_PriceTheta_list = ST_data_1[:, 3]
# 股权市场价值列表
ST_Plus_E_list = ST_Plus_Clpr_list * ST_Plus_Lsttrdshr_list + (
        ST_Plus_Trdshr_list - ST_Plus_Lsttrdshr_list) * ST_Plus_NAPS_list
# Short-term Debt (in 10,000)
ST_Plus_SD_list = ST_Plus_data_2[:, 0]
# print(ST_SD_list)
# Long-term Debt (in 10,000)
ST_Plus_LD_list = ST_Plus_data_2[:, 1]
# print(ST_LD_list)
# 负债合计
ST_Plus_D_list = ST_Plus_data_2[:, 2]

# 非ST公司数据
# 收盘价_Clpr
NST_Clpr_list = NST_data_1[:, 0]
# 流通股_Trdshr
NST_Trdshr_list = NST_data_1[:, 1]
# 已上市流通股_Lsttrdshr
NST_Lsttrdshr_list = NST_data_1[:, 2]
# 每股净资产
NST_NAPS_list = NST_data_1[:, 5]

# 无风险列表
NST_r_list = NST_data_1[:, 4]
# 股权收益率波动率列表
NST_PriceTheta_list = NST_data_1[:, 3]
# 股权市场价值列表
NST_E_list = NST_Clpr_list * NST_Lsttrdshr_list + (NST_Trdshr_list - NST_Lsttrdshr_list) * NST_NAPS_list
# Short-term Debt (in 10,000)
NST_SD_list = NST_data_2[:, 0]
# print(ST_SD_list)
# Long-term Debt (in 10,000)
NST_LD_list = NST_data_2[:, 1]
# 负债合计
NST_D_list = NST_data_2[:, 2]


def KMVOptSearch(E, D, r, T, EquityTheta):
    """
       KMV模型计算
    """
    EtoD = float(E) / float(D)

    # print("EtoD：", EtoD)
    # print()
    # global daishu
    # daishu = 1

    def KMVfun(x, EtoD, r, T, EquityTheta):
        """
            解d1,d2方程组
        """
        EtoD, r, T, EquityTheta = float(EtoD), float(r), float(T), float(EquityTheta)
        d1 = (np.log(x[0] * EtoD) + (r + 0.5 * x[1] ** 2) * T) / (x[1] * np.sqrt(T))
        d2 = d1 - x[1] * np.sqrt(T)
        # global daishu
        # print("第{}次计算：".format(daishu))
        # daishu += 1
        # print(x[0], x[1])
        # print(x[0] * stats.norm.cdf(d1, 0.0, 1.0) - np.exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0) / EtoD - 1,
        #       stats.norm.cdf(d1, 0.0, 1.0) * x[0] * x[1] - EquityTheta
        #       )
        # print()
        return [
            x[0] * stats.norm.cdf(d1, 0.0, 1.0) - np.exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0) / EtoD - 1,
            stats.norm.cdf(d1, 0.0, 1.0) * x[0] * x[1] - EquityTheta
        ]

    # print()
    VaThetaX = fsolve(KMVfun, [1, 0.1], args=(EtoD, r, T, EquityTheta))
    Va = VaThetaX[0] * E
    AssetTheta = VaThetaX[1]
    return Va, AssetTheta


def calculate_DD(i, counter, r_list, SD_list, LD_list, a, b, D_list, PriceTheta_list, E_list, IS_ST_bool):
    """
        计算DD距离并返归m，n值
    """
    m = 0
    n = 0
    for time in range(counter):
        # 利用KMV模型计算
        r = r_list[time]  # 无风险利率
        T = 1  # 期权到期时间，以年为单位
        SD = SD_list[time]  # SD为长期负债 （单位：元）
        LD = LD_list[time]  # LD短期负债 （单位：元）
        DP = a[i] * SD + b[i] * LD  # 违约触发点DPT # DP=（TD-LD） +0.5LD （7） 其中：TD为总负债；LD为长期负债。 （单位：元）
        D = D_list[time]  # D=TD-LD+0.5∗LD  债券的面值是D （单位：元）
        PriceTheta = PriceTheta_list[time]  # 交易波动率，可以是日月年的数据 （单位：%）
        EquityTheta = PriceTheta  # 股权收益率的波动率σE
        E = E_list[time]  # 公司股票的市场价值 （单位：元）
        Va, AssetTheta = KMVOptSearch(E, D, r, T, EquityTheta)  # 资产价值VA（单位：元）资产收益率波动率σA
        AssetTheta = AssetTheta  # 资产收益率波动率σA （单位：%）
        DD = (Va - DP) / (Va * AssetTheta)  # 违约距离DD
        EDF = stats.norm.cdf(-DD)  # 预期违约概率EDF（单位：%）
        # print("ST_DD = {}".format(DD))
        if IS_ST_bool == 1:
            if DD >= 0:
                m += 1
        elif IS_ST_bool == 0:
            if DD < 0:
                n += 1

    return [m, n]


# Attention: you had better put the aim-function in another file.
def aim(Phen, LegV):  # define the aim function
    """
        定义目标函数，输出表现型和可行矩阵
    """
    global counter
    counter += 1
    print("第{}次迭代开始：".format(counter))

    print(Phen)

    for i in range(Phen.shape[0]):
        # print("N = ", N)
        n = 0
        m = 0
        a = Phen[:, 0]
        # print("a:", a)
        b = Phen[:, 1]
        # print("b:", b)

        # 计算ST公司的m，n值
        mn_list = calculate_DD(i, count, ST_r_list, ST_SD_list, ST_LD_list, a, b, ST_D_list, ST_PriceTheta_list,
                               ST_E_list,
                               1)

        m += mn_list[0]
        n += mn_list[1]

        # 计算*ST公司的m，n值
        mn_list = calculate_DD(i, count, ST_Plus_r_list, ST_Plus_SD_list, ST_Plus_LD_list, a, b, ST_Plus_D_list,
                               ST_Plus_PriceTheta_list, ST_Plus_E_list, 1)
        m += mn_list[0]
        n += mn_list[1]

        # 计算非ST公司的m，n值
        mn_list = calculate_DD(i, num_com, NST_r_list, NST_SD_list, NST_LD_list, a, b, NST_D_list, NST_PriceTheta_list,
                               NST_E_list, 0)
        m += mn_list[0]
        n += mn_list[1]

        print("i:", i + 1)
        print("a = {}, b = {}".format(a[i], b[i]))
        print("m = {}, n = {}".format(m, n))
        # print(1 - ((m + n) / N))
        Mis_rate = 1 - ((m + n) / N)
        print("Mis_rate：", Mis_rate)
        if i == 0:
            diff = np.array([[Mis_rate]])
        else:
            diff = np.append(diff, [[Mis_rate]], axis=0)
        print("diff.shape:", diff.shape)
        print()

    print("[diff, LegV][diff, LegV][diff, LegV][diff, LegV]:", [diff, LegV])
    return [diff, LegV]


def main():
    """
        主函数：主要设定遗传算法变量和调用
    """

    # AIM_M = __import__('aimfuc')
    AIM_M = __import__('main2')  # get the handle of aim-function
    # variables setting

    # 变量设置
    x1 = [0, 10]  # 自变量1的范围
    x2 = [0, 10]  # 自变量2的范围
    b1 = [0, 1]  # 自变量1是否包含下界
    b2 = [0, 1]  # 自变量2是否包含上界

    ranges = np.vstack([x1, x2]).T  # 生成自变量的范围矩阵
    # print(ranges)
    borders = np.vstack([b1, b2]).T  # 生成自变量的边界矩阵
    # print(borders)
    # 生成区域描述器
    # help(ga.crtfld)
    FieldDR = ga.crtfld(ranges, borders)  # create FieldDR
    # print(FieldDR)

    # help(ga.sga_new_real_templet)  # 查看 ga.sga_new_real_templet函数的调用方法以及参数设置
    # call the GEA algorithm template
    [pop_trace, var_trace, times] = ga.sga_new_real_templet(AIM_M, 'aim', None, None, FieldDR, problem='R', maxormin=1,
                                                            MAXGEN=200, NIND=40, SUBPOP=1, GGAP=0.9,
                                                            selectStyle='etour',
                                                            recombinStyle='xovdprs', recopt=0.9, pm=None,
                                                            distribute=True, drawing=1)

    # output results
    for num in var_trace[np.argmin(pop_trace[:, 1]), :]:
        print(chr(int(num)), end='')


# 若调用该文件则调用main()主函数
if __name__ == "__main__":
    main()
